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We classify the ground states and topological defects of a rotating two-component condensate when 
varying several parameters: the intracomponent coupling strengths, the intercomponent coupling 
strength and the particle numbers. No restriction is placed on the masses or trapping frequencies 
of the individual components. We present numerical phase diagrams which show the boundaries 
between the regions of coexistence, spatial separation and symmetry breaking. Defects such as 
triangular coreless vortex lattices, square coreless vortex lattices and giant skyrmions are classified. 
Various aspects of the phase diagrams are analytically justified thanks to a non-linear cr-model that 
describes the condensate in terms of the total density and a pseudo-spin representation. 

PACS numbers: 



I. Introduction 

Bose Einstein condensates (BECs) provide an excel- 
lent environment to study experimentally and theoreti- 
cally a rich variety of macroscopic quantum phenomena. 
In a rotating single component condensate, topological 
defects to the order parameter often manifest themselves 
as vortices that correspond to a zero of the order pa- 
rameter with a circulation of the phase. When they get 
numerous, these vortices arrange themselves on a trian- 
gular lattice. In fact, vortices were first observed in two- 
component BEC's [Ij. Since then two-component BECs 
and the topological excitations within have been experi- 
mentally realised in a number of configurations: a single 
isotope that is in two different hyperfine spin states [l|-ul ^ 
two different isotopes of the same atom [8[ or isotopes of 
two different atoms 0-ll|. 

When a two-component condensate is under rotation, 
topological defects of both order parameters are cre- 
ated which lead to more exotic defects such as singly or 
multiply quantised skyrmions 0, [H - 16 1 . Analogy with 
the Skyrme model from particle physics is often invoked 
to represent the defects [isj . The singly quantised 
skyrmions contain a vortex in one component which has 
the effect of creating a corresponding density peak in the 
other component. These singly quantised skyrmions are 
often referred to as coreless vortices. Once numerous, 
these vortices and peaks arrange themselves in (coreless 
vortex) lattices, that can be either triangular or square. 
Other defects such as vortex sheets or giant skyrmions 
can also be observed l3"21|. The aim of this paper is 
to classify the type of defects according to the different 
parameters of the problem restricting ourselves to two- 
dimensions. While this paper will only concern itself with 
magnetically trapped two-component BEC's, there is ac- 



tive research into spinor condensates (for a recent review 
see [H). 

In the mean-field regime, the two-component Bose- 
Einstein condensate at zero temperature is described 
in terms of two wave functions (order parameters), 
^^i and \E'2, respectively representing component- 1 and 
component-2. The two-component condensate is placed 
into rotation about the z-axis with il = Clz where Cl is 
the rotation frequency assumed to apply equally to both 
components. The Gross-Pitaevskii (CP) energy func- 
tional of the rotating two-component, two-dimensional 
BEC is then given by 



(1) 



' 2/^, and Vfc(r) = mfea;^r^/2 are the har- 
monic trapping potentials with trapping frequencies Wk, 
centered at the origin. Here mi and m2 are the masses 
of the bosons in component- 1 and component-2 respec- 
tively, and h is Planck's constant. The angular momen- 
tum is in the z-axis and is defined &s ~ i[z-{rx p)] for 
linear momentum p. The energy functional ([T|) contains 
three interaction constants: Uk representing the internal 
interactions in component k, and U12 representing the 
interactions between the two components. 

The time-dependent coupled Gross-Pitaevskii (GP) 
equations are obtained from a variational procedure, 
ihd'^k/dt = SE/S'^l, on the energy functional ([T]). 
Let oj — (wi + Ld2)/2 be the average of the trap- 
ping frequencies of the two components and introduce 
the reduced mass TO12 such that — itii^ + m^^. 
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The GP energy and the coupled GP equations can be 
non-dimensionahsed by choosing Cj^^, Huj and rg = 
^/h/ (2?7ii2w) as units of time, energy and length respec- 
tively. On defining the non-dimensional intracomponent 
coupling parameters gk = 2Uk'mi2/fi^ and the intercom- 
ponent coupling parameter 512 = 2Ui2mi2/h? (= 321), 
the dimensionless coupled GP equations read 



.dipk 



(V - lAk) ipk + -. 

rrik 4toi2 



for = ri/cj and where 



A, 



1 

2 mi2 



Jl X r, 



(2) 
(3) 



for J7 = (0,0,0) and r = (x, j/,0). The ground state 
of the energy ([T|) or the GP equations ^ is determined 
by preserving the normalisation condition which in this 
paper is either taken to be 



(fr ■ 



Nk 



where Nk is the total particle number of the kih- 
component, or 



(l^l|' + l^2|' 



(5) 



The following sections will only consider repulsive in- 
teractions so that gi, g2 and 1712 are always non-negative. 
In order to separate the regimes of interest, a non- 
dimensional parameter combining these gk and gi2 will 
appear: 



ri2 



1 



gf2 
.9152' 



(6) 



Furthermore to simplify matters, we introduce ratio pa- 
rameters "q and ^ such that 



mi — r/m2 and uJi ~ ^lu2 with {?7,^} > 0. 



(7) 



The effective trapping potentials for each component 
coming from equation ^ are then, respectively, 



Vf{r) ={ij + 1) 
_(77+l) 



(^ + 1)^ 
1 



(8a) 
(8b) 



It follows that the limiting rotational frequency for each 
component is O'l™ = 2C/(^ + 1) and O^™ = 2/(^ + 1), 
so that it is necessary to consider a rotational frequency 
such that 0<n< O'™ = min{0'/"\ O^""}. 



Experimentally, it is often the case that the C/fc, nik, ujk 
and Nk are of the same order, so that much of the theo- 
retical and numerical analysis concerning two-component 
condensates has assumed equality of these parameters. 
In the case where the intrac omp onent coupling strengths 
are equal, Kasamatsu et al. |l9l.l20j produced a numerical 
phase diagram in terms of the rotation and the intercom- 
ponent coupling strengths: they found phase separation 
regions with either vortex sheets or droplet behavior and 
regions of coexistence of the components with coreless 
vortices, arranged in triangular or square lattices. 

In this paper, we present a classification of the ground 
states and various types of topological defect when the 
intracomponent coupling strengths are distinct. Depend- 
ing on the magnitude of the various parameters the com- 
ponents can either coexist, spatially separate or exhibit 
symmetry breaking. A richer phenomenology of topolog- 
ical defects is then found. 

Much of the analysis carried out to investigate the 
ground states and topological defects will use a a non- 
linear sigma model. It has been introduced previously in 
the literature [l^, 24 1 for = ^ = 1 but we will generalise 
this to the cases rj and ^ different from 1. This involves 



(4) writing the total density as 



PT = |V'l|' + ^K'2|' 



(9) 



, X2]^ is in- 



A normalised complex- valued spinor x = [xii 
troduced from which the wave functions are decomposed 
as -01 = y/pTXi, 1P2 \/pt/vX2 where |xip -I- |x2p = 1- 
Then we define the spin density S — X^^X where <t are 
the Pauli matrices. This gives the components of S as 



X =XlX2 +X2XI, 



Si, 



i{x*iX2 - X2X1), 



S.=\xi\'-\X2\ 



(10a) 
(10b) 
(10c) 



with jSp = 1 everywhere. We will write the GP energy 
([1]) in terms of px and S that will allow us to derive 
information on the ground state of the system. 

The paper is organised as follows. The different re- 
gions of the n — ri2 phase diagrams are outlined in Sect. 
II with a detailed analysis of the range of ground states 
and topological defects. Then the non-linear cr-model is 
developed in Sect. Ill to analyse the ground states in 
terms of the total density. Lastly Sect. IV takes the sec- 
ond normalisation condition ([5]) and presents an example 
O — ri2 phase diagram. 



II. The ri2 - n Phase Diagram 

A. Numerical Parameter Sets 

In this paper, three different configurations are consid- 
ered, two of which relate directly to experimental config- 
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urations. Firstly, we analyze a ^^Rb-^''Rb mixture with 
one isotope in spin state \F = 2, nif = 1) and the other 
in state |1, —1). The complex relative motions between 
these two isotope components of rubidium were consid- 
ered experimentally by Hall et al Here, the masses 
are equal (77 = 1) and the transverse trapping poten- 
tials are equal, wi =102 = 1). The scattering lengths 
for each component are oi = 53.35A and 02 = 56.65A. 
The intracomponent coupling strengths Uk for a two- 
dimensional model are defined as Uk = VSnll^ak/liTikazk] 
(k = 1,2), and the intercomponent coupling strength is 
defined as U12 = ^/2TTh^ai2/[mi2az], where mi2 is the 
reduced mass, given by = "^-T^ + ''^2^ ^ '^k, 0-12 are 
the s-wave (radial) scattering lengths, azk is the char- 
acteristic size of the condensate in the z direction, and 
a-z = (flzi + az2)/2. In the experiment the computa- 
tions of approximate values of Qzi and az2 can be made 
using a Thomas-Fermi approximation in the z direction 
as explained in [25|, chapter 17. This relies on the as- 



sumption that NkUk/cizk are large. Then, it is reasonable 
to describe the experiments through a two-dimensional 
model. For our simulations, we choose gi = 0.0078 and 
52 = 0.0083, which are consistent with the experimen- 
tal data and we set Ni ^ N2 ^ 10^. We denote this 
parameter set as 'Experimental Set 1' (ESI). 

Secondly, we tackle a ''^K-^^Rb mixture with both iso- 
topes in spin state |2, 2), as was considered by Modugno 
et al 0. Component-1 is identified with the '^^K iso- 
tope and component-2 with the *^Rb isotope. In the ex- 
periment, the masses are different (77 ~ 0.48), but since 
77^^ = 1, then miiof = 77i2w| so that both components 
experience the same trapping potential. Here the scat- 
tering lengths for each component are ai ~ 31.75A and 
a2 = 52.39A and we choose intracomponent coupling 
strengths as gi = 0.0067 and 52 = 0.0063. This set is 
denoted ES2. 

Lastly, in order to make an analogy with previous the- 
oretical studies, we consider the set (ES3) which contains 
a mixture chosen such that all the parameter groups are 
equal, i.e. here gi = 52 = 0.0078 and Ni = N2 = 10^ 
with equality of the 771^ and uJk {r] = ^ = 1). 

In each case, the features of the ground states will be 
explained as gi2 and ft are varied. While the value of the 
gk are nominally fixed, the value of the intercomponent 
coupling strength, gi2, can be altered by Feshbach reso- 
nance (see for instance 2^2^), which allows for an ex- 
tensive experimental range in ri2, given by ([6]) (keeping 
ri2 < 1)- Simulations are thus performed on the coupled 
GP equations ^ in imaginary time. In general, it is dif- 
ficult to find the true minimizing energy state. But the 
use of various initial data converging to the same final 
state allows us to state that the true ground state will be 
of the same pattern as the one we exhibit. 

We would like to note that while we have chosen these 
particular parameters, we have conducted extensive nu- 
merical simulations over a range of different parameters 



sets and find the sets presented here illustrate well the 
possible ground states. 



B. Classification of the Regions of the Phase 
Diagram 

The ground states can be classified according to 

1. the symmetry properties 

2. the properties of coexistence of the condensates or 
spatial separation. 

When there is no rotation (fJ = 0), as illustrated in 
Fig. [U the geometry of the ground state can either be 

• two disks (coexistence of the components and sym- 
metry preserving state) 

• a disk and an annulus (symmetry preserving with 
spatial separation of the components), which de- 
pends strongly on the fact that the gk, ruk, uJk are 
not equal 

• droplets (symmetry breaking and spatial separa- 
tion). 

Under rotation, the different ground states can be clas- 
sified according to the parameters ri2 and VI. When the 
condensates are two disks, or a disk and annulus, de- 
fects may break some symmetry of the system as Q, is in- 
creased. But the wave function retains some non trivial 
rotational symmetry. We will not refer to this as symme- 
try breaking since the bulk condensate keeps some sym- 
metry. We will use the terminology symmetry breaking 
when the bulk completely breaks the symmetry of the 
system and is a droplet or has vortex sheets. We find 
that there are four distinct regions, determined by the 
geometry of the ground state. 

Region 1: Both components are disk shaped, with no 
spatial separation. Above some critical velocity fi, core- 
less vortices appear: a vortex in one component which 
has the effect of creating a corresponding density peak 
in the other component. They arrange themselves either 
on a triangular or on a square lattice. Figure [2] provides 
a typical example of the density profiles. 

Region 2: Vortex sheets. Under the effect of strong ro- 
tation, the components spatially separate and completely 
break the symmetry of the system. Nevertheless, they 
are approximately disk shaped with similar radii. Many 
vortices are nucleated that arrange themselves into rows 
that can have various patterns: either they can be striped 
or bent and are often disconnected from each other. A 
vortex sheet in one component corresponds to a region 
of macroscopic density in the other component. These 
features can be seen in the density plots of Fig. EJc). 
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FIG. 1: Density (divided by 10*) plots along y — for 
component- 1 (dashed lines), component-2 (dotted lines) and 
the total density (solid lines) in which the components (a) 
coexist both being disks (ri2 = 0.1), (b) spatially separate to 
be a disk and an annulus (ri2 = 0) and (c) have symmetry 
broken to be two droplets (ri2 — —0.3). The parameters are 
gi = 0.0078, 52 = 0.0083 and iVi = iV2 = lO'^ with 77 = ^ = 1 
(set ESI). The angular velocity of rotation is f2 = 0. Distance 
is measured in units of ro and density in units of r,^^. 



Region 3: Spatial separation preserving some symme- 
try. Here one component is a disk while the other com- 
ponent is an annulus and the disk fits within the annulus 
with a boundary layer region in which both components 
have microscopically small density as illustrated in Fig- 
urelSl^b) anddljb). Under rotation, the topological defects 
can either be coreless vortex lattices in the disk and (but 
not always) corresponding isolated density peaks in the 
annulus, and/or a giant skyrmion at the boundary inter- 
face between the two components, as will be described 
later. 

Region 4- Rotating droplets. The components spa- 
tially separate and have broken symmetry such that the 
centres of each condensate arc different and each compo- 
nent contains a single patch of density as illustrated in 
Figures [3lja) and|4lja). In the rotating droplets, vortices 
can appear provided the rotation is greater than some 
critical value. These features can be seen in the density 
plots of Fig. E 

The f2 — phase diagrams corresponding to the ex- 
perimental parameter sets introduced above (sets ESI, 



FIG. 2: (Color online) A series of density plots for component- 
1 (left column) and component-2 (right column) in which the 
components coexist and are both disks. The parameters are 
gi = 0.0078, 32 = 0.0083 and iVi = iV2 = 10^ and with 
77 = C = 1 (set ESI) and ri2 = 0.5 (which gives gw = 0.0057). 
The angular velocity of rotation is Q and it takes the values 
(a) 0.25, (b) 0.5, and (c) 0.75. Vortices in one component 
have a corresponding density peak in the other component 
(coreless vortices). In (c) the coreless vortex lattice is square. 
At these parameters the components are in region 1 of the 
phase diagram of Fig. [6l Distance is measured in units of ro 
and density in units of r^'^. 



ES2 and ESS) are shown in Fig. 's El [7] and [H The last 
one is of similar type as the one reported by Kasamatsu 
et al ^ (there gi = g2 = 4000 and TVi = TVs = 1/2). 

New features can be observed in Fig [6] and [7l such 
as isolated density peaks that eventually vanish as 
Ti2 is made more negative and the multiply quantised 
skyrmions at the interface between the disk component 
and the annular component in region 3. When {i],^} ^ 1 
(Fig. [7]), no region 2 is found to exist. This absence is 
easily explained by two factors: the onset of region 3 
for (large) positive values of and the lack of vor- 
tices nucleated in component-1 (and to some extent in 
component-2). 

We will now analyze in more detail some of the features 
of regions 1, 3 and 4 of the phase diagrams. Vortex sheets 
(region 2) have been analysed (mainly in the case of equal 
gk, rrik and Wfe) in 
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FIG. 3: (Color online) A series of density plots for component- 
1 (left column) and component-2 (right column) in which the 
components have spatially separated. The parameters are 
gi = 0.0078, 52 = 0.0083 and Ni = N2 = 10^ and with 
r; = 5 = 1 (set ESI) and = —0.3 (which gives gi2 — 
0.0092). The angular velocity of rotation is and it takes the 
values (a) 0.1, (b) 0.5, and (c) 0.9. In (a) the components are 
rotating droplets (region 4 of Fig. [6]), in (b) the components 
have spatially separated (but keep some symmetry) and have 
isolated density peaks (region 3 of Fig. [6]) and in (c) there 
are vortex sheets present (region 2 of Fig. [6|. Distance is 
measured in units of ro and density in units of r^'^. 



C. Analysis of the Features of the Phase Diagrams 

1. Square lattices 

It is a specific property of two components to stabi- 
lize square lattices. A requirement for the existence of 
square lattices is the nucleation of many vortices in both 
components, something not permitted in ES2 (Fig. [7|), 
where only a small number of vortices are ever nucleated 
in component 1. 

Square lattices generally occur at high rotational ve- 
locities and examples have been observed experimentally 
and numerically [3, [H [H. Mueller and Ho (si 



and later Kegeli and Oktel [33[ have analysed the tran- 
sition from triangular to square lattices as is var- 
ied in two-component condensates when gi ^ 172 and 
is such that the condensate is in the lowest Landau 
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FIG. 4: (Color online) A series of density plots for component- 
1 (left column) and component-2 (right column). The pa- 
rameters are gi = 0.0078, 52 = 0.0083 and iVi = 7V2 = 10^ 
and with ri = ^ = 1 (set ESI) and ri2 = —1.3 (which gives 
gi2 = 0.0122). The angular velocity of rotation is SI and it 
takes the values (a) and (b) 0.9. In (a) the components are 
rotating droplets (region 4 of Fig. [S]) while in (b) the com- 
ponents have spatially separated (but keep some symmetry) 
and there are no isolated density peaks (region 3 of Fig. 
Distance is measured in units of ro and density in units of 

_-2 



level (LLL). Providing fl is large, according to [32 



the value of a = \/l — T12 determines whether the vor- 
tex lattices are triangular (0 < a < 0.172), in transi- 
tion to becoming square (0.172 < a < 0.373) or square 
(a > 0.373). Therefore, assuming that, in the example 
parameters of gi = 0.0078, 32 = 0.0083 and 77 = ^ = 1 
(set ESI), the gk arc sufficiently close enough, the square 
lattice should first appear at a = 0.37, or cquivalently 
ri2 = 0.86. Comparing this to the numerically ob- 
tained values, where the square lattice first appears for 
0.8 < ri2 < 0.9 in Fig. |6l the agreement seernsgood. It 
will be interesting to see how the analysis of 3^, l33[ has 
to be modified when gi becomes distant from 32- 

Nevertheless, the square lattices are present for lower 
rotational velocities at which the LLL is not necessarily 
valid. In a forthcoming paper, using the nonlinear sigma 
model presented below, we hope to derive a vortex energy 
in terms of the positions of vortices. The ground state 
of this energy indeed leads to square lattices in some 
range of parameters (see also Section III.E). It turns out 
that this vortex energy is similar to that of Barnett et al 

Mm. 
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FIG. 5: (Color online) A density plot for component-1 (left 
column) and component-2 (right column) showing examples 
of vortex nucleation in rotating droplets. The parameters are 
gi = 0.0078, 32 = 0.0083 and iVi = iV2 = lO'^ and with 
■q = ^ = l (set ESI) and = -9 and Q, = 0.5. At these pa- 
rameters the components are in region 4 of the phase diagram 
of Fig. [6l Distance is measured in units of ro and density in 
units of Tq '^. 



2. Rotating droplets 

In the parameter range of region 4, symmetry break- 
ing with spatial separation occurs. When the conden- 
sate is not under rotation, the ground states of the two 
components will be 'half ball'-like structures, as can be 
seen in Fig. IDJa) (if gi = 32 and -q — — 1, then the 
two components are exactly half- balls). The difference 
between the intcrcomponcnt strengths introduces a cur- 
vature to the inner boundary of each component with the 
component that corresponds to larger having positive 
curvature and the other component having negative cur- 
vature. This curvature depends on both and 17. If 
is held constant, then the curvature increases as 17 
increases. When ri2 goes to —00, the droplets approach 
half-balls. 

Vortex nucleation is also seen in region 4; see Fig. [S] 
In this figure there are four vortices in component-1 and 
three vortices in component-2. The number of vortices 
in each component will increase as 17 is increased but 
this increase will also increase the curvature of the inner 
boundaries of the components, thus preventing the vor- 
tices aligning themselves into vortex sheets. Examples 
of the rotating droplet ground states can be seen in Fig. 
m^a), and further examples can be found in [22.. i36i i39i] . 



3. Spatial separation preserving some symmetry: disk plus 
annulus 
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FIG. 6: (Color online) ^ — T\2 phase diagrams for parameters 
gi = 0.0078, 52 = 0.0083, iVi = iV2 = 10^ with n = i = I 
(set ESI) with normalisation taken over the individual com- 
ponents (Eq. (Q). (a) Numerical simulations where trian- 
gles (squares) indicate that the vortex lattice in both com- 
ponents is triangular (square) and diamonds that no vortices 
have been nucleated. Filled triangles, squares and diamonds 
are where the two components are disk-shaped and coexist. 
The empty triangles and empty diamonds are where the two 
components have spatially separated: one component is a disk 
and the other an annulus with a giant skyrmion at the bound- 
ary layer; those triangles with a dot in the centre represent 
the appearance of coreless vortices inside the annulus. The 
crosses 'x' are where the two components have broken sym- 
metry and are vortex sheets and the crosses ' + ' are rotating 
droplets, (b) A schematic representation of the numerical 
simulations. The solid lines indicate the boundary between 
different identified regions (determined by the geometry of 
the ground state) and the dashed lines the boundary between 
triangular and square lattices in region 1 and the boundary 
between peaks and no peaks in region 3. The boundary be- 
tween region 1 and the others can be calculated analytically 
by Eq. ((39)) to give ri2 = 0.008. The unit of rotation is Cj. 



Region 3 is defined by the ri2 and U in which the 
condensates have spatially separated components but 
still possess some symmetry about the origin: the disk 
component-1 sits securely within the annular component- 
2 (see Fig.'s[31Jb) and|Hb)). There is a boundary layer 
evidenced where the outer edge of the disk overlaps the 



inner edge of the annulus. 

Let us describe the onset of region 3 from region 4, 
captured in Fig. [3] (sub- image (a) to (b) or (c) to (b)). 
For a particular ri2, as 17 is increased, then the curvature 
increases to such an extent that the components develop 
constant non-zero curvature, and are identified as a disk 
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Region 3: 
Spatial separation (disk 
and annuius) 




Giani skyrmions, possibility of a single ring 
of vortices and/or density peaks in 
component-2. / 


/, Region 1 : 
/ Coexistence 
(2 disks) 




Triangular coreless vortex 
iattices or rings of vortices. 


Region 4: Rotating / 
droplets \ 


^ ^ No vortices 



-0.6 -0.4 -0.2 0.2 

ri2 



FIG. 7: (Color online) n — Ti2 pliase diagrams for parameters 
gi = 0.0067, g2 = 0.0063, Ni = N2 = lO" with = 0.48 and 
77^^ = 1 (set ES2) with normalisation taken over the individ- 
ual components (Eq. Q). (a) Numerical simulations where 
triangles indicate that the vortex lattice in component-2 is 
triangular, circles that component-2 contains rings of vortices 
and diamonds where no vortices have been nucleated. Filled 
triangles, circles and diamonds are where the two components 
are disk-shaped and coexist. The empty circles and empty 
diamonds are where the two components have spatially sep- 
arated: one component is a disk and the other an annuius 
with a giant skyrmion at the boundary layer; those circles 
with a dot in the centre represent the appearance of coreless 
vortices inside the annuius. The crosses are where the two 
components have broken symmetry and are rotating droplets, 
(b) A schematic representation of the numerical simulations. 
The solid lines indicate the boundary between different identi- 
fied regions (determined by the geometry of the ground state). 
The boundary between region 1 and region 3 is also calculated 
analytically by Eq. (|39|l (dotted line). For these parameters 
there is no region 2. The unit of rotation is uj. 



and an annuius. Conversely, if f2 is held constant and 
is pushed to more negative values, then the curvature 
decreases. 

Under rotation, defects can be observed: coreless vor- 
tices and giant skyrmions. Coreless vortices sit inside the 
disk while giant skyrmions are observed in the boundary 
layer. 



a. 0.6 




Square coreless 
vortex lattices 



Region 1 : 
Coexistence 
(2 disks) 



Triangular coreless 
vortex iattices 



ri2 



FIG. 8: A schematic ^l~Ti2 phase diagram for gi — g2 = g = 
0.0078, iVi = iV2 = 10^ where 77 = ^ = 1 (set ESS) and with 
the normalisation taken over the individual components (Eq. 
(lU). The solid lines indicate the boundary between different 
identified regions (determined by the geometry of the ground 
state) and the dashed lines the boundary between triangular 
and square lattices. For these parameters there is no region 
3. The unit of rotation is Cj. 



4- Coreless Vortiees in the disk plus annuius case 

In the geometry of disk plus annuius, the vortices in 
the inside disk have corresponding isolated density peaks 
inside the annuius, hence in the microscopic density re- 
gions, as illustrated in Fig. |3|b). Thus any vortex lat- 
tice structure produced in component- 1 is replicated by 
a peak lattice structure in component-2, i.e. there is the 
continuing presence of coreless vortex lattices in the spa- 
tially separated condensates (ioj . 

As the number of vortices in the disk increases with 
increasing angular velocity, the number of density peaks 
inside the central hole of the annuius likewise increases. 
This has the effect that the central hole of the annuius 
can be masked at high angular velocities by the multiple 
occurrence of the density peaks. A recent analytical un- 
derstanding of the interaction between vortices and peaks 
has been provided by 4l| . It may be extended to the case 
of the disk and annuius, for which the average density in 
the central hole is very small for one component. 

Pushing to lower values has the effect of reduc- 
ing the size of the boundary layer between the disk and 
annular components, but also the isolated density peaks 
disappear; see Fig. IH^b). A recent work has analysed, 
from an energy perspective, the preference for the ground 
state to contain or not contain density peaks [42| . 

Figure [9] plots the maximum peak density of the den- 
sity peaks in component-2 as a function of for the 
parameters of Fig. [H] and with ft = 0.5 and = 0.75. 
The disjoint region when il = 0.75 for -0.4 < < 0.1 
on Fig. ini is due to the absence of density peaks in 
component-2 as a result of the appearance of vortex 
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FIG. 9: Plots of maximum peak density in component-2 
against ri2 for Q = 0.5 (dashed line) and Q = 0.75 (solid 
line) when gi = 0.0078, 52 = 0.0083, Ni = N2 = 10^ and 
with rj — = 1 (set ESI) . The disjoint region when Q, = 0.75 
for —0.4 < < 0.1 is due to the absence of isolated density 
peaks in component-2 as a result of the appearance of vortex 
sheets. Density is measured in units of r^^. 



sheets (region 4). We see from Fig. [SI that the maxi- 
mum peak density occurs when ^ 0.01 which is the 
value at which the components begin to spatially sepa- 
rate (the transition between region 1 and region 3) in set 
ESI. For and that take values outside of region 
1, the maximum density of the peaks decreases linearly 
as T12 increases. The maximum density approaches mi- 
croscopically small values for — 1.1 when ft = 0.5 
and for ^ —1 when f2 = 0.75. An example of the 
transition can be seen in Fig.'s[3Jb) and|H[b). 



5. Giant skyrmion 

A boundary layer between the overlap of component- 1 
and component-2 is present for all values of ri2 in region 
3 but reduces in width as ri2 becomes more negative (in- 
deed the boundary layer disappears only as ri2 — > —00). 
There are additional topological defects that occur in the 
boundary layer that are not discernible on a traditional 
density plot. These topological defects' presence can be 
observed in a phase profile, however a better visualisa- 
tion is to use the pseudospin representation (P))- ((TU)) . One 
can plot the functions Sx, Sy and/or Sz which reveal the 
presence of all the topological defects - the coreless vor- 
tices (singly quantised skyrmions) that have already been 
visualised on the plots, and the new defect, a multiply 
quantised skyrmion. 

The distinct nature of the two types of topological de- 
fect can be clearly seen in Fig. [TOl a.b) where an Sx plot 
and an {Sx,Sy) vectorial plot are shown respectively. A 
density plot of each component for these same parame- 
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FIG. 10: (Color online) Plots of (a) Sx and (b) a vectorial plot 
of iSx,Sy) for Fia = -0.3 and = 0.5 when gi = 0.0078, 
ff2 = 0.0083, iVi = iV2 = 10^ and with r? = ? = 1 (set ESI). 
The ring of skyrmions, at r = rs, traces the boundary layer 
and the coreless vortices exist for r < Vs- Blow-ups (c) and 
(d) of (b) highlight the nature of the two types of topological 
defect; (c) a coreless vortex at (1, —1) and (d) a section of the 
multiply quantised skyrmion. The respective density plots 
for this case are shown in Fig. [S^b). Distance is measured in 
units of vq. 



ters is seen in Fig. Elb). The coreless vortices evident 
in Fig. ^h) are again clearly evident in Fig. [TUf a) and 
(b). A blow-up of the region close to one of (the three) 
coreless vortices in the {Sx,Sy) vectorial plot, centred 
at (1,-1) and exhibiting circular disgyration is shown 
in Fig. fTUTc). The texture of S can also exhibit cross- 
and radial- disgyration Conversely, the multiply 

quantised skyrmion, not present in the density plots of 
Fig. ^h), can clearly be seen in both Fig. [TPT a) and 
(b) . The multiply quantised skyrmion forms a ring along 
the boundary layer. A blow-up around (0, —3.5) for the 
{Sx, Sy) vectorial plot again shows the multiply quantised 
nature of this defect. A multiply quantised skyrmion was 
evidenced by [2lj who termed it a giant skyrmion, a ter- 
minology retained in this paper. Increasing the rotation 
results in an increase in both the number of coreless vor- 
tices and in the multiplicity of the giant skyrmion. Ref 
2lj gave a relationship, m = n + 1, between the multi- 
plicity of the giant skyrmion n and the total circulation 
m in the central hole of component-2. However the nu- 
merical simulations here suggest that a more appropriate 
correlation is provided hy m = n + I where / is the num- 
ber of isolated coreless vortices. As an example, consider 
Fig. [TO] where = 0.5 and n = 8, Z = 3 and m = 11, 
satisfying m — n + l. It should be noted that the relation- 



ship m = n + I still holds as the isolated coreless density 
peaks disappear from the density profile of coniponent-2; 
this is because the phase imprints due to the vortices in 
component-1, which constitute I, are still present. 



III. A Non-Linear a-Model 

The ground states of the rotating two-component con- 
densate can be recovered from a 'non-linear tT-model' for- 
mulation of the energy functional in terms of the total 
density pr and spin vector S (Eq.'s ([9|)-(fT0|)). Wc will 
analyze in detail the possible regimes in this section and 
show how a Thomas Fermi approximation can hold in the 
case Ti2 > or can be generalized and provide relevant 
information in the case T12 < 0. 



A. Energy Functional Representation 

We write the energy functional of the two wave func- 
tions (Eq. dD), non-dimensionalised) as i?[V'i,V'2] = 
Eke + EpE + Ei where 



E^ 



1 



KE 



V-^(r/ + l)f^xr)^/'i 



E^ 



PE 



r] + l 
f ?/ 



V - — (?7 + l)n X r ) ^^2 



(11a) 



(lib) 



Ei = j ^ffilV-il' + \92M + gi2M\i^2 



with 



2 (i + e)2 
1 (1 + ^) 



1 



]2 



277(1+^2 8^ 



(He) 

(12a) 
(12b) 



It is assumed that the g^, m^j, w^, and iV^ are distinct 
so that a weighted total density can be defined as © 
with ■01 = y^prXi and 'ip2 = V Pt/vX2- The spin den- 
sity vector S, which satisfies = 1 everywhere, has 
components given by Eq.'s ([T0| . We then have 



I^i|' = ^Pt(1 + 5,), \^,\^ = ±p^(1-S,_). (13) 

We introduce the phases 9i and 62 defined by xi = 
\Xi\exp{i6i) and X2 = IX2I exp(i02). The energy func- 
tional is expressed in terms of 4 variables: the total den- 
sity pt, the component Sz and the angles di, 02- We see 



that Eke = Ep^ + Esz + Eg^,g^ where 
1 



EpT 

Es. 
Eei,e2 



—AV^Y d'r. 



iv + 1) 

1 PT (V5,)^ 
4(77 + 1) (1-^2) 

1 PT 



Sr. 



(14a) 
(14b) 



{i + Sz) ( ve*! - -{l + ri)n X r 



2(1 + 77) 



+ {l-Sz)[ye2-—{l + ri)nxr 



(fr. 



(14c) 



The other terms of the energy straightforwardly become 



E, 



PE 



[(ji + nh) + (ji - ]2h)Sz] r^pT dli^&) 



Pt 



Ei^ / ^ (^0 + ^1^- + ^2^,') d\ (15b) 



with 



Co 



Cl 



1 
1 
1 



(jy^gi + .92 + 277312), 
{Tfgi - 92), 



C2 =^(^7 91 + 92- 277^12). 



(16a) 
(16b) 
(16c) 



Thus the complete energy is 



E = 



SI) 



Pt 



2(1 + 77) 



(l + 5,)( V0i--(l + 7,)nxr 



(17) 



+ {l^Sz)[Ve2-—{l + 'n)^y.r 
+ [(ji + ^2/??) + (ji - 32N)Sz\ r^PT 



Pt 



Co + ciS'^ + C2SJ d r. 



Energy ([TT]) is subject to the constraints that can be 
rewritten as 



Pt(1-5'^)/2 d^r^N2T], (18a) 
PTil + Sz)/2 d^r^Ni. (18b) 



10 



or equivalently 



B. A Thomas- Fermi approximation 



PT (fr =Ni + N2ri, 
PtS^ (fr =Ni - N2r], 



(19a) 
(19b) 



The minimization of the energy under the constraints 
(Uni) yields two coupled equations with two Lagrange mul- 
tiphers, p and A. We wiU write them when the phase is 
constant (the gradient of the phases 9i and 02 can be 
ignored) and the effect of rotation is contained in an ef- 
fective centrifugal dilation of the trapping potential in- 
cluded in ji,j2- 



fi + XSz 



1 



1 



(V5, 



(1+77) 4(1 + 77) (1-52) 

+ [in + nh) + (ji - ]2h)Sz\ 

+ Pt(co + cyS-, -f C2SI), 



(20) 



and 



1 



4(1 + 7;) (1-^2) 2(1 + 7,) (I-52) 

+ (ji - 32hy + (ci + 2c25,) pt/2. 

(21) 

As pointed out in [l^ , in the general case where 77 and 
^ are different from 1, we have seen from the expression 
Eg-^ g^ that the effective velocity in each component is 
different. Nevertheless, in the case when 77 = ^ = 1, it 
is possible to define an effective velocity shared by both 
components, 



"off ^ 



ve Rs, 

~ 2(1 - 52) 



(22) 



where Q ^ di 
the identity 



and R = SyVS^ 



S^VSy. We note 



(23) 



where {VSf = {VS^^ + (VSy)^ + [VS.f. Expansion 
of the square in Eg-^_g^, substituting in the 7;off and using 
the identity from above reduces the energy to the simple 
form found in 



^ {v,s -nxrf + ^r\l-n^)pT (24) 

+ ~^ {^0+ ClSz + C2SI 



where cq, ci and C2 arc equal to cq, ci and C2 with 77 set 
equal to unity. 



The typical Thomas-Fermi (TF) approximation to 
Eq.'s (P(I| - (|2ip is to assume that derivatives in px and 
Sz are negligible in front of the other terms. If we apply 
the TF approximation we then get 



p + A^, = [(ji + J2/77) + Ui - j2/v)Sz] 

+ Pt(co + CiSz + C2Sz), 



and 



A = (ii - 32/vy + ^ (ci + 2025*^) PT- 
The TF energy is then 

Etv = I [{ji +j2/v) + (ji - j2hl)Sz] r^PT 



(25) 



(26) 



Pt 



Co + ciSz + C2S^) ~ PPT ~ XprSz dr. 

(27) 



It is important to point out that the reduction of this 
quadratic form in px and pxSz yields 



^TF = / ^ (ptSz + ^PT + ^{ji - 32/^1^ - ^ 

2 \ 2C2 C2 C2 



2 Ko 



4C2 



1 



Pt 



(co - C?/4C2) 

ci 



1 

'2(coC2-c2/4) 
+ constant terms. 



{p-{h+]2h)r') 
2 



(A - (ji + .72/??)r") 



d^r 



(28) 



The existence of a minimum for this quadratic form 
implies that C2 > and cq — cf/(4c2) > 0. Since 
cq — c\/ (4C2) = gig2^i2/ (47;2c2), this implies in particular 
that ri2 > 0. 

Therefore, if ri2 > 0, (which implies C2 > 0), a 
Thomas Fermi approximation can be performed as has 
previously been considered, generally taking an approx- 
imation on the individual component wave functions t/^i 
and '02 [H, [36, Usij. Nevertheless, as we will see be- 
low, if C2 < and ri2 < 0, then (cq - ^) > and a TF 
approximation can still be performed on p^, provided we 
keep gradient terms in Sz- 

Multiplying (^5)) by Sz and subtracting (^5)) leads to 
P = (Ji + n/vV + PT [^Sz + Co) • (29) 
Simplification by ptSz with (j26p yields, when t/ii x ^j2 7^ 
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181) then give 



and 



where 



Pt 



03 + a^r 

9192^12 



03 + a4r2 



ai = 477 (Aco - /iCi/2), 

02 = »75i^2 - 52^11, 

03 = 477^ (/iC2 - Aci/2), 

04 = -{r]gih2 + g2hi), 



and where we define 



Jk 



.912 



(30) 



(31) 



(32a) 
(32b) 
(32e) 
(32d) 



(33) 



When only one component is present (tpi x ■02 = 0), this 
simphfies to 



U[, + X-2,nr^] if 02=0 
^ UMM-A)-2j2r2] if ^1=0 



since ~ +1 when -02 and S2 = —1 when 0i = 0. 

To begin, we note from Eq. ([31]) that = (oi — 
aj,Sz)/{aASz - 02), so that 



Pt 



aiQi — 0203 
gi92^i2{aASz - 0,2) 



(35) 



Since the intracomponent couphng strengths are chosen 
to be distinct, the support of each component wiU not 
necessarily be equal. In order to lead the computations, 
we have to assume a geometry for the components: two 
disks, a disk and annulus or two half-balls (droplets) and 
a sign for ri2. 



C. > 



1. Both components are disks 



27r 



(1 + S.) 
Pt rar 







(0104 - 0203) 
4:9192^12 



-1 



(1 + S,) 
^ {ttiSz - 02)^ 



(36a) 



where sq = Sz\r=o ~ ai/as and 



VN2 
2ti 



(1 - , 
PT 7: rdr 

^ 



(aia4 - 0203)^ 

4:9192^12 



R2 

PT rdr 

Ri 

{aiSz - a2Y 



dS, 



92 

4r]j2 



Pt dpT, 



with 



Pd =- 



as + a4r^ 



r=Ri 



9192^12 



0.20.3 — oiai 



(36b) 



(37) 



since Sz 



9192^12(02 + 04) ' 

-1 at r = 



Completion of these integrals and noting that {02 
04, 02 - 04} = {-2g2hi,2rigih2} gives 



2 _ SNigiglTi2hi 



(38a) 



7r(so + 1)2 

which must necessarily be positive, i.e. hiTi2 > 0, and 

7Vi(a2 - soa4)2 =2r]'^ N29ig2i2hiTi2{l + sq)^ 
+ 27Vi(l + so)7?.72,9iri2X 

[7751/12(1 + So) + 2,92/11 (so - 1)] 

(38b) 



where in Eq. (j38bp the expression for 03 from Eq. ([38a| 
has been substituted. This equation can always be solved 
in terms of so when hiTi2 is positive since the discrim- 
inant is equal to SNxrf g29ij2hiTi2{Nigi2 + N292)- We 
find that 



1+so 



2iVi52 



^1.92 - r]Ni9i2 + \^'2Niri'^gij2Ti2{Nigi2 + A^252)/^i 



We start when both components are disks. Without 
loss of generality, we can assume that the outer bound- 
ary of component-2 (at r = R2) is larger than that of 
component-1 (at r = Ri). Therefore, Sz ~ —1 and ([34]) 
holds in the annulus, while psp holds in the coexisting 
region, which is the inside disk. The integrals from Eq.'s 



A similar calculation can be completed if the outer 
boundary of component-2 is larger than that of 
component-1. This then gives that h2Ti2 > (although 
expressions (|38|) change slightly). 

If an annulus develops in component 3 — fc, it means 
that So = 1 (so = —1) for k = 1 (k = 2). Inputting this 
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choice into Eq. (j38bp gives 



that 



5l2 



2{Niji + N2]2) 



Nkgkjs-k 




Nljl + N2j2 



4jV3-fcgig2j3-> 

Nljl + N2j2 



1/2 (39) 



as the critical gi2 at which an annulus forms in 
component-{3 — k}. In the case of Fig. [71 the curve 
has been plotted in dashed lines and is close to the 
numerical curve. Notice that if = 1 (equal trapping 
frequencies for both components), then Eq. ([M]) becomes 
independent of J7, as in the phase diagrams of Fig. [51 
where it yields = 0.008. 

Let us point out that before the transition to the disk 
plus annulus takes place, there is a subregion of region 1, 
where there are 2 disks, but in one component the wave 
function has a local minimum at the origin. For instance, 
in the case of Fig. |6l it corresponds to Fi changing sign. 

As a conclusion, in order for the ground state to be 
composed of two disks, assuming that componcnt-fc is the 
component with smaller support, we need that hk > 0, 
> and gi2 < gi2- These three conditions can be 
summarised as 



5-12 < mm 



■^'^ gs^k, ^^gk,gi2, V5i.92 ) if/ii, /i2 > 0, 



33-k 



3k 



(40a) 



j3-fc 

jk 



gk < gi2 < min 



Jk 
js-k " 



g3-k,gi2, V51.92 ) if/i3-fe<0 

(40b) 



2. A disk and an annular component 

In this case, we can assume a disk in component-1 
and an annulus in component- 2; there are three regions: 
an inner disk where only component-1 is present and 
P4p holds, an outer annulus where only component- 2 is 
present and (j34p holds and an inner annulus where both 
components coexist and psp holds. 



In order to use the TF approximation when one com- 
ponent is annular, computations similar to P8ap - (j38bp 
lead to hih2 < and gi2 > gi2- This can be summarized 
(for an annulus in component- {3 — k}) as 

f -^3-/0 - \ ^ ^ ■ ( I 

max — — gk,g\2 < 3i2 < mm gs-fe, V5i52 

\ 3k ) \33-k 

(41) 

with /i3_fc < 0. Equation (HIT) also places the restriction 



gk <g3-k 



3k 



j3-k 

=.93-feAfc, 



(42) 



where = {jk/j3-k)'^- Note that A1A2 = 1 such that an 
annulus develops in component-2 (-1) if (72 > (<) gi^2- 



3. Orders of Intracomponent Strengths and Special Cases 

The effect that changing the order of the intracom- 
ponent strengths and particle numbers has on gi2, and 
thus on the phase diagrams, is now investigated. There 
are two cases to consider, depending on the relative or- 
ders of the particle numbers. Drawing aid from experi- 
mental values, it is always expected that minjA^i, 3> 
max{f/i, (72}. Throughout it will be assumed that g2 > 
A2gi, so the annulus develops in componcnt-2 and that 
ji and j2 are of order unity. 

In the first case when the particle number of 
component-1 is much greater than the particle number of 
component-2 (iVi N2), it follows that gi2 ~ gi32/ ji- 
The boundary between region 1 and region 3 is then di- 
rectly dependent on the value of the ratio J23i/[ji52]- 
Notice that if g2 ^ gi, T12 evaluated at 312 = 512 tends 
to unity and as such an annulus would always be present 
in component-2, whatever the value of 312. 

Conversely in the second case when the particle num- 
ber of component-1 is much smaller than the particle 
number of componcnt-2 (A^i ^ A^2), it follows that 
512 ~ V5152 which implies that the annulus will only 
develop near ri2 = 0. 

We can also look at some special cases - there are four 
that can be considered: 

Case (i). Akgs-k = gk- When A2 (equivalently Ai) is 
such that this equality is made, the two components are 
both disks and no annulus develops. 

Case (a). ?/ = ^ = 1. Then A^ = 1 and R2 ^ i?i 
A^252ri ^ A^i^iFi when there are two disks and the an- 
nulus develops in the component which has the larger 
interaction strength. 

Case (in). '7 = C — 1 ^-nd 5i = 52 = 5- When the 
intracomponent coupling strengths are equal, gi2 — g 
and there are always two disks with R2 ^ i?i iV2 ^ 

Ni. 

Case (iv). rj = ^ = 1, gi = g2 ^ g and Ni = N2 = N. 
When the particle numbers and intracomponent coupling 
strengths are equal, §12 = g and there are always two 
disks with R2 ~ Ri- A detailed analysis of this case, 
explored numerically in a phase diagram for all ri2 and 
analytically in the TF limit, was considered by [20j . 
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4- Justification of the Thomas-Fermi approximation 

The gradient terms in (pn)) -(pT |) can be neglected if d^i, 
the characteristic length of variation of p-r and is much 
smaller than dc, the characteristic size of the condensates. 
We have that l/d^j is of order of /x and A, which are of 
order ^/N^gk, \/Nkg3-k- Hence dhi is bounded above by 
the maximum of {Nkgk)~^^'^, {Nkg3-k)~^^^ ■ From the 
expression of 03, the characteristic size of the condensate 
is of order the minimum of {gkNkri2y/\ Therefore, 
the Thomas Fermi approximation holds if NkgkV^i2, 
Nkg3-kV^i2 are large. This requires the usual Thomas 
Fermi criterion that Nkgk, Nkga-k are large, but breaks 
down if is too small. 

For = 0, the equations lead to spatial separation: 
either the radii of the disks tend to in the case of two 
disks or the outer radius of the inner disk tends to the in- 
ner radius of the annulus in the case of disk plus annulus 
so that ipi X 'ip2 = everywhere. Another analysis has to 
be carried out to understand the region of coexistence, 
which is of small size and has strong gradients. 



(a) (b) 
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FIG. 11: Total density profiles (divided by 10*) obtained 
numerically (solid lines) and analytically (dashed lines) for 
two spatial separation cases with gi — 0.0078, 32 ~ 0.0083, 
Ni ^ N2 = 10^ ry = C = 1 (set ESI) and T12 = -10 and 
n = 0: (a) annulus plus disk, analytical estimate coming 
from (|43|l and (b) droplet, analytical estimate coming from 
(|49|) . The inset in (a) shows the discontinuity of density at 
Ts = 3.63. Distance is measured in units of ro and density in 
units of r^^. 



D. < 0, beyond the TF approximation 

For negative T12, and if Nkgk are large, the Thomas 
Fermi approximation can be extended provided some 
model takes into account the small region where the con- 
densates coexist. 



Indeed, if we go back to ([28|) . and have both C2 < and 
< 0, then the coefficient in front of the second square 
is positive and the optimal situation is to have the square 
equal to 0, which leads to the inverted parabola (pO| . On 
the other hand, the coefficient in front of the first square 
is negative, and the ground state involves derivatives in 
Sz to compensate it. Under the assumption that the 
boundary layer where Sz varies is small, we are going to 
derive a TF model with jump for px- We will analyze it 
for the different geometries (disk plus annulus, droplets 
and vortex sheets) and show that it provides information 
consistent with the numerics. 



1. Disk Plus Annulus 



Assuming that the boundary layer is present only at 
some r = Ts, then S'^ = +1 in the region in which 
component-2 is zero (r < rj) and 5z = — 1 in the re- 
gion in which component- 1 is zero (r > r+). The transi- 
tion from = 4-1 to 5*2 = — 1 is not smooth, therefore 
creating the jump in density. Therefore, we are lead to 



minimise the integral 

I = J^ 2jir2pT + d\ 



2 , 52 2 72 

2— r pT + tTzPt d r 



(43) 



with respect to Tj. Here Br^ is a ball of radius and 
BR\r, is a torus with outer boundary at r = i? and inner 
boundary at r = Tg. Thus 



27rr,, 



o ■ 2 I .91 2 
2jir PT + —Pt 



2nr. 



o-?2 2 
2—r PT 



92 2 



(44) 



= 0, 



which implies that 
1 



p =— [p,i - 2jir^] r e Br^ 
91 

=— [m2 - 2j2r^] r e BR\rs- 
92 



(45a) 
(45b) 



Then using the normalisation conditions J p d^r = Ni 
and J p^d^r = r\N2 we get an outer radius 



(46a) 
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and chemical potentials 
Nigi 



Ai2 =2j2 




(46b) 
(46c) 



It remains to find r^. But, from (|44|) . it follows that 



Jl 



Ni 



N2J2 

■^92 



(47) 



The above leaves a quartic equation in r^, the solution 
of which can be found numerically. 

We can then complete the energy calculation to give 



91 N? 
27rr2 



-iNdi+2N2j2ys 



6gi " 3 \ TT 



1/2 



(48) 



A check plot of the density profile with the parameters 
gi = 0.0078, 52 = 0.0083 and iVi = TVa = 10^ where 77 = 
^ = 1 agrees well with this model where the boundary is 
calculated to be at rg = 3.62 (sec Fig[TTJa)). 



2. Droplets 

This formalism can also be extended to study the 
droplet case. As before, assuming a thin boundary layer, 
one can allow for a jump in density. We thus need to 
minimise the Thomas-Fermi energy for two droplets de- 
scribed by Bi and B2 that exist in the regions < 9 < a 
and a < 6 < 2tt respectively. The energy is then 



Ri 

(•R2 



gi 2 



rdr 



(49) 



V 



The expressions for the density in each domain (j45p and 
the normalisation conditions allow completion of this in- 
tegral to give 



4^/2 
~3~ 



N^gui 



1/2 



Nig2j2 
(27r - a) 



1/2- 



(50) 



It remains to find the optimum a. This is achieved 
through the condition dl /da = which gives 



(A^fsiJi)^/^ {Nlg2j2)'/' 



v3/2 



a 



{2tt - a)3/2 



(51) 



27r 



(1 + Nifjm 



where we have set N = N1/N2, g = gi/g2 and j = 
ill 32- As expected, equality of the iV^, gu and setting 
rj = ^ = 1 (=> jl ~ J2) gives a = tt and the condensate is 
then composed of two half-balls. Otherwise, a curvature 
is present. A check plot of the density profile with the 
parameters gi = 0.0078, 32 = 0.0083 and iVi = 7V2 = 10^ 
where 77 = ^ = 1 is given in Fig lllf b) . 

Finally we can note the energy for the droplets is 



1 



30F 

\NlgmY" 



, \ -1/2 

{NiNi{giglnf2)"'f" 



(52) 



This energy can be compared to the energy of the disk 
plus annulus (|48p to determine which is the optimum 
geometry. Indeed, in the numerical cases studied before, 
the droplets are preferred states for small Vl. 



3. Regions of Vortex Sheets 

In the case of vortex sheets, we can assume that the 
global profile of the total density is TF-like, obeying Eq. 
([50]) in the bulk of the condensate. By working with the 
total density, we do not require any information on the 
vortex sheets themselves (and consequently Sz)- 

We thus take the form of pr from Eq. ([5(7| from 
which we note that the outer boundary at r = i? satisfies 
R — \/ —azja^ and that completion of the normalisation 
condition (|19ap gives 



03 = 



2(iVi +777V2)5ig2ri2a4 



(53) 



This expression evidently requires a4ri2 < 0. We how- 
ever expect the vortex sheets to be present only in the 
ri2 < domain, thus we can be more specific on the 
condition, namely that 04 > 0, or 



512 > 



??gij2 + g2ji 

{Vjl+j2) 



(54) 



We point out that this critical number corresponds to 
ri2 = in the cases studied above. With as above 
and a4 given by the parameters of the system, the density 
profile is then fully accessible. 



E. Analysis of defects 

The advantage of the Thomas Fermi analysis of px in 
the nonlinear sigma model is that it allows for analysis of 
defects as a perturbation calculation when Nkgk, Nkga-k 
and Nkgi2 are large, in the spirit of what has been done 



for a single condensate in [44|, l52H54f or for two conden- 
sates in 



We do not need to analyze specifically the 
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peaks because they are taken into account in the Sz for- 
mulation. This has been attacked in a more complicated 
fashion in the case of a single coreless vortex in [19| . We 
will develop our ideas in a later work but let us point out 
that the starting point is the energy (jl7p and if we call pi 
and Qi the location of the vortex in each component, the 
main terms coming from the vortex contributions lead 
to: 



Each 



of each component provid- 
energy term proportional to 
+ 77)) where dhi is the healing 



vortex core 
ing a kinetic 
Pt(0) log(iw/(2(l 
length, the characteristic size of the vortex core, 
here of order 1 / \/Nkgk- 

• the rotational term providing a term in 
— cr2(pT(0))^, where c is a numerical constant. The 
balance of these two terms allows us to compute 
the critical velocity for the nucleation of the first 
vortex. 

• the kinetic energy then yielding a term in — log \pi~ 
Pj\ - log|% - 

• the rotation term yielding a term in r2(|pip + 

• the interaction term providing a term p^S^ for the 
perturbation of Sz close to a coreless vortex. The 
ansatz can be made in several ways leading to an 
interaction term in e~'^'~**' . 

This leads to a point energy of the type 

E«(IP'l' + l'?^l'+^e-l^''-9'l')-^(log|K-p,|+log|g.-g,| 

where a and b are related to the parameters of the prob- 
lem. The ground state of such a point energy leads to 
a square lattice for a sufficient number of points and for 
some range of a and b. 



IV. Phase Diagram Under Conserved Total Particle 
Number 

In the experiment of Hall et al Q a single compo- 
nent EEC of *^Rb in the |1,-1) state was initially cre- 
ated, before a transfer of any desired fraction of the 
atoms from this |1,— 1) state to the |2,1) state created the 
two-component BEC. Thus the ratio of particle numbers 
-^1 /N2 is controllable experimentally, with the constraint 
that Ni + N2 is constant (in the case of the experiment 
of m , A^i -f- A^2 = 5 X 10^). As such, experimentally, it 
is possible to keep the individual particle numbers con- 
stant (as in the normalisation condition (|4])) or to keep 
the total particle number A^i -I- A^2 constant, allowing Ni 
and N2 to vary (as in ([5])). We produce an O — phase 
diagram for the parameters gi = 0.003, 52 = 0.006 and 
77 = ^ = 1 using the normalisation condition given by 



a. 
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FIG. 12: (Color online) Q, — T\2 phase diagram for parameters 
gi = 0.003 and 32 = 0.006 with 77 = ^ = 1 using the nor- 
malisation condition given by ([5]). (a) Numerical simulations 
where triangles indicate that the lattice in both components 
is triangular, squares that the lattice in both components is 
square and diamonds where no vortices have been nucleated. 
Filled triangles, squares and diamonds are where the two com- 
ponents are disk-shaped and coexist, empty triangles, squares 
and diamonds are where only component-1 exists; those with 
a dot in the centre represent the appearance of coreless vor- 
tices in component-2 and those without a dot in the centre 
represent the complete disappearance of component-2. (b) A 
schematic representation of the numerical simulations. The 
solid lines indicate the boundary between different identified 
regions (determined by the geometry of the ground state) and 
the dashed lines the boundary between triangular and square 
lattices. The unit of rotation is Cj. 



(0 with Ni+ N2 = 2.1 X 10^. Note that if gi = g2, 
mi = TO2 and oji = a;2 (?? = C = l)i then the phase 
diagram would be identical to that of Fig. [8] (i.e. the 
normalisation condition in this case would not be impor- 
tant). The f2 — ri2 phase diagram is presented in Fig. 
[T2] There are three distinct regions (determined by the 
geometry of the ground state): 

Region 1. In the first region, both components are disk- 
shaped. As before, the coreless vortices can either form 
a triangular or a square lattice depending on the values 
of ri2 and ri. Figure [T3] shows this case for ri2 = 0.7 
and VL = 0.65 (where a triangular lattice is present) and 
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(a) (a) 




FIG. 13: (Color online) A series of density plots for 
component- 1 (left column) and component-2 (right column). 
The parameters are gi = 0.003 and g2 = 0.006 with rj = — 1 
and ri2 — 0.7 (which gives gi2 — 0.0023) and normalisation 
taken over the total density (Eq. The angular velocity 

of rotation is SI and it takes the values (a) 0.65 and (b) 0.9. 
There is a triangular lattice in (a) and a square lattice in (b). 
At these parameters the components are in region 1 of the 
phase diagram Fig. 1121 Distance is measured in units of ro 
and density in units of t-q^. 



il = 0.9 (where a square lattice is present). In region 1, 
both components have the same radii. 

Region 2. In the second region, only component- 
1 exists except for isolated density peaks that exist in 
component-2. These isolated density peaks occur at the 
same location as the vortices do in component- 1, and are 
thus identical to the isolated coreless vortices described 
in detail in Sect. IV. Figure [HI shows this case ~ 0.3 
with fl = 0.5 and n = 0.9. 

Region 3. In the third region, only component- 1 exists 
(the density peaks in component-2 that were present in 
region 2 are no longer present). Furthermore, only tri- 
angular vortex lattices are observed in this effective one 
component condensate (in component-1). 

Computations similar to Section III hold, except that 
now we have to set A = 0. This leads to ai/a^ = 
— Ci/(2c2). We see that this ratio (which is Sz{Q)), 
reaches 1 or -1 when ci = ±2c2. In our numerical 
case, this leads to f 12 = 0.5. We see clearly 3 regimes: 
ri2 > f 12, where the condensates are 2 disks, ri2 < 0, 
which is phase separation, in which case S'^ = 1 is the 
preferred state and < ri2 < f 12, in which case the TF 
approximation leads to a computation of pT with core- 
less vortex lattices and variations in Sz which improve 
the energy and lead to this intermediate state, still to be 



FIG. 14: (Color online) A series of density plots for 
component-1 (left column) and component-2 (right column). 
The parameters are gi = 0.003 and g2 ~ 0.006 with 77 = ,^ = 1 
and — 0.3 (which gives gi2 = 0.0035) and normalisation 
taken over the total density (Eq. ([5}). The angular velocity 
of rotation is Q, and it takes the values (a) 0.5 and (b) 0.9. 
At these parameters the components are in region 2 of the 
phase diagram Fig. 1121 Distance is measured in units of ro 
and density in units of r^^. 

studied in more detail. 



V. Conclusion 

We have presented phase diagrams of rotating two 
component condensates in terms of the angular velocity 

and a nondimcnsionnalized parameter related to the 
couphng strengths ri2 = 1 — S'i2/(s'i52)- We have ana- 
lyzed the various ground states and topological defects 
and have found four sets characterized by the symme- 
try preserving/symmetry breaking, coexistence or spatial 
separation of the components. When the geometry of the 
ground states is either two disks (coexistence of compo- 
nents, region 1) or a disk and an annulus (spatial sepa- 
ration keeping some symmetry, region 3), the topological 
defects are coreless vortex lattices (with possible stabi- 
lization of the square lattice) or giant skyrmions at the 
boundary interface between the two components. In the 
complete symmetry breaking case, we have found vortex 
sheets and droplets. The difference of masses or coupling 
strengths between the components can induce very dif- 
ferent patterns. 

We have introduced an energy pT|l related to the total 
density and a pseudo spin vector. The minimization in a 
generalized Thomas Fermi approximation provides a lot 
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of information on the ground states for general masses, 
trapping frequencies and coupling strengths. Some parts 
of the phase diagrams can be justified rigorously, both 
in the case > 0, which had been studied before, but 
also in the case < with generalized models. This 
formulation of the energy should bring in the future more 
information on the defects. 
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